Gao et at. BMC Genomics 2013, 14:738 
http://www.bionnedcentral.conn/1471 -21 64/1 4/738 



Genomics 



RESEARCH ARTICLE Open Access 



Do transgenesis and marker-assisted backcross 
breeding produce substantially equivalent 
plants? - A comparative study of transgenic and 
backcross rice carrying bacterial blight resistant 
gene Xa21 

Lifen Gao^'^^ Yinghao Cao^'^^ Zhihui Xia\ Guanghuai Jiang\ Guozhen Liu^ Weixiong Zhang^'^" 
and Wenxue Zhai^" 



Abstract 

Background: The potential impact of genetically nnodified (GM) plants on human health has attracted much 
attention worldwide, and the issue remains controversial. This is in sharp contrast to the broad acceptance of plants 
produced by breeding through Marker Assisted Backcrossing (MAB). 

Results: Focusing on transcriptome variation and perturbation to signaling pathways, we assessed the molecular 
and biological aspects of substantial equivalence, a general principle for food safety endorsed by the Food and 
Agricultural Organization and the World Health Organization, between a transgenic crop and a plant from MAB 
breeding. We compared a transgenic rice line (DXT) and a MAB rice line (DXB), both of which contain the gene 
Xq21 providing resistance to bacterial leaf blight. By using Next-Generation sequencing data of DXT, DXB and their 
parental line (D62B), we compared the transcriptome variation of DXT and DXB. Remarkably, DXT had 43% fewer 
differentially expressed genes (DEGs) than DXB. The genes exclusively expressed in DXT and in DXB have pathogen 
and stress defense functions. Functional categories of DEGs in DXT were comparable to that in DXB, and seven of 
the eleven pathways significantly affected by transgenesis were also perturbed by MAB breeding. 

Conclusions: These results indicated that the transgenic rice and rice from MAB breeding are substantial equivalent 
at the transcriptome level, and paved a way for further study of transgenic rice, e.g., understanding the chemical 
and nutritional properties of the DEGs identified in the current study. 
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Background 

The primary objective of breeding in agriculture is to de- 
velop plants of desired genotypes or traits, such as high 
yields and resistance to adverse environmental impact. 
Marker-assisted backcrossing (MAB) and transgenesis 
(aka genetic modification or GM) are two widely adopted 
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plant breeding techniques. As a conventional technique, 
MAB breeding has been used to develop new crop culti- 
vars of, e.g., baley, maize and rice [1-4]. The basis of MAB 
breeding [5] is to transfer a specific allele at the target 
locus in a donor line to a recipient line while selecting 
against donor introgression across the rest of the genome. 
In most cases, the recipient line used for backcrossing has 
a large number of favorable attributes but is deficient in a 
few characteristics. Since MAB breeding to certain degree 
mimics or replicates natural selection, novel cultivars pro- 
duced through MAB breeding have been regarded as 
genetically safe. However, MAB breeding is laborious. 
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requiring several backcrossings and a large number of 
individual plant screenings, typically on the order of thou- 
sands. It typically takes a great deal of luck to produce a 
product of a desirable trait. 

Transgenesis is an effective means for improving crop 
genetic makeup for deriving favorable traits. Breeding by 
transgenesis has several advantages over breeding by 
MAB. Transgenesis is a direct means for introducing a 
gene or genotype to a genome in order to produce a tar- 
get trait. As a result, the potential to produce plants with 
favorable traits increases dramatically. For example, by 
transfering a gene that encodes a plastidial adenyalte 
kinase into a potato variety, the transgenic potato displays 
wild-type growth and developmental phenotype but also 
has an increased yield and starch concentration [6]. By 
introducing the hybrid cellulase gene cel-hybl into a 
spring barley variety through Agrobacterium-mediated 
transformation, the selected marker-free transgenic barley 
produces a high level of cellulase (l,4-|3-glucanase) in de- 
veloping grains, suggesting that the transgenic barley has 
the potential for producing a large quantity of cellulase for 
commercial use [7]. The nutritional value of Golden Rice 
is improved with increased pro-vitamin A content by 
introducing genes encoding phytoene synthase (psy) in 
combination with the Erwinia uredovora carotene de- 
saturase {crti) into rice [8]. Through an Agrobacterium- 
mediated genetic transformation system, Xa21, a rice 
bacterial blight resistance gene, has been introduced into 
five Chinese rice varieties and as a result, the transgenic 
rice plants exhibit a high resistance to bacterial blight [9]. 

While transgenesis offers immense opportunities to cur- 
tail the severe threat of food shortage the expanding world 
population is facing, there are considerable public con- 
cerns over the use of transgenesis for crop improvement. 
Indeed, it remains extremely controversial whether or not 
transgenic crops have an adverse impact on human health. 
At the center of this controversy is the issue of whether or 
not insertion of a transgene into the host plant genome or 
manipulation of an allele in the host genome may affect 
the expression of other genes and ultimately lead to unin- 
tended phenotypes. Unfortunately, it is technically chal- 
lenging to address this issue because accurate prediction 
of phenotypes based on genotypic variation and/or gene 
expression alteration remains a research topic. 

Due to the difficulty, efforts of evaluating the safety of 
transgenic crops has been geared toward assessing the 
substantial equivalence between transgenic plants and 
wildtype or conventionally bred plants, like plants from 
MAB breeding [10]. Substantial equivalence has been in- 
troduced as a standard by the Organization for Economic 
Cooperation and Development (OECD) and has been 
endorsed by the Food and Agriculture Organization of the 
United Nations/World Health Organization (FAO/WHO) 
[11]. However, the standard is based on comparative 



analysis and offers only a general principle. No specific 
molecular, biological, chemical or nutritional basis has 
been established to precisely specify the degree of substan- 
tial equivalence [12,13]. It thus leaves widely open the 
study of various aspects of equivalence, ranging from mo- 
lecular, biological, and chemical to nutritional equivalence, 
between a transgenic plant and a wildtype or plant pro- 
duced by MAB breeding. Nevertheless, it has been agreed 
that to be considered substantially equivalent, the charac- 
teristics of a transgenic plant must be within a natural 
range of variation [14], a guideline we follow in our study. 

Rice is an essential staple crop for the world population 
and a model plant for basic and applied research. Rice 
bacterial leaf blight (BLB), caused by bacteria Xanthomonas 
oryzae pv. oryzae {Xoo), is one of the most devastating rice 
diseases throughout the world. Utilization of BLB resistant 
genes in breeding is the most effective and economical 
strategy for controlling BLB [15]. Xa21, the first-cloned 
BLB resistant gene from Oryza longistaminata, has re- 
ceived much attention because of its broad spectrum of re- 
sistance to BLB [16]. The gene has also been widely used in 
BLB resistance breeding through both the transgenic and 
MAB strategies [17-21]. Therefore, transgenic rice and 
MAB rice carrying Xa21 offer an excellent opportunity to 
assess the possible substantial equivalence of transgenic 
and MAB rice as well as rice in natural environments. 

In order to pave the way for future studies of the safety of 
transgenic crops, we focused on the molecular and biological 
aspects of substantial equivalence of transgenic rice. We 
adopted a systems -biology perspective and examined the 
transcriptome variation of transgenic rice. Specifically, we in- 
corporated Xa21 into the three-line maintainer line D62B 
through transgenesis and MAB breeding. Adopting Next 
Generation (Next-Gen) sequencing, we profiled the tran- 
scriptomes of four rice plants: the Xa21 transgenic line 
(named as DXT), the Xa21 MAB breeding line (named as 
DXB), the untransformed recipient D62B, and another rice 
variety MH86 (restorer line). We then analyzed transcrip- 
tome variation of the two rice plants carrying Xa21 in refer- 
ence to that of D62B and transcriptome change between 
D62B and MH86. This transcriptome analysis was further 
enhanced by a pathway analysis to understand the pathways 
that might be disturbed in the two rice plants carrying Xf^2i. 

Results 

A system for comparative study of transgenic 
and MAB rice 

In order to compare transgenic and MAB plants, transgenic 
and MAB rice plants carrying Xa21 using the parental line 
D62B were constructed. To generate the transgenic rice, 
Xa21 was introduced into D62B through Agrobacterium- 
mediated transgenesis. The transgenic rice plants were se- 
lected from the Tq, Ti and T2 generations by molecular and 
resistance analysis. The homozygous, single copy, and 
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marker-free transgenic line, or DXT for short, was obtained 
from the T2 generation. To confirm the results. Southern 
hybridization analysis was performed using the restriction 
endonuclease Pvull to digest genomic DNA from the trans- 
genic and parental plants. The details for developing the 
transgenic rice was described in our early report [22]. The 
sample used for analysis was the T9 generation of DXT with 
stable agronomic traits. 

Xa21 was introgressed into the parental line D62B to 
produce the MAB breeding line using IRBB21 as the 
donor. IRBB21 was bred by transfering Xa21 into IR24 
through backcrossing [23]. Six backcrossing generations 
were made because it is usually necessary to take a mini- 
mum of six backcrossing generations in order to recover 
the phenotype of recurrent parent lines and eliminate 
donor chromosome fragments linked to the target gene 
[5]. A backcrossed line with homozygous Xa21 and simi- 
lar phenotype with the recipient D62B was obtained in 
BC6F2 generation and named DXB. 

In order to facilitate direct in-field screening and molecular 
analysis of transgenic and MAB plants that showed consist- 
ent agronomic traits similar to that of their parental line, the 
transgenic line (DXT), the MAB line (DXB) and their paren- 
tal line (D62B) were grown in the same fields in the breeding 
process. D62B can thus serve as an ideal control for the com- 
parison of DXT and DXB. In order to introduce a reference 
to natural variation, another rice varieties MH86, an indica 
restorer line in the three-line breeding system was also in- 
cluded in the profiling experiments. Since rice carrying Xa21 
confers robust resistance to most strains of Xoo at adult 
stages, the RNA samples were extracted from adult leaves of 
the four rice lines for transcriptome profiling. 

The transgenic rice and MAB rice were 
phenotypically similar 

The morphological characteristics of DXB, DXT, D62B, 
and MH86 were examined in the rice fields. DXT and 
DXB were morphologically similar to their parental line 
under visual inspection. The major agronomic traits of 
DXT and DXB, listed in Table 1, were also scored in the 
fields. The results showed that the main agronomic traits 
of DXT and DXB were consistent with or comparable to 
that of D62B, except with respect to the focal trait of 
BLB resistance. Both DXT and DXB were highly resistant 



Table 1 Major agronomic traits of four rice plants studied 



Sample 


Lesion length (cm) 


Plant height (cm)^ 


Tillers 


D62B 


10.0±2.7 


88.4±4.3 


11±3 


DXT 


0.4±0.3* 


86.3±2.0 


10±1 


DXB 


1 .2±0.4* 


85.5±3.4 


12±2 


MH86 


11.3±1.9 


11 3.1 ±2.8* 


10±2 



to nine Philippine races of Xoo but D62B was susecptible 
to the all Xoo races (Figure 1). Highly phenotypic simi- 
larity among DXT, DXB, and D62B suggested that both 
transgenic and MAB breeding strategies had very little 
impact on the morphological characteristics and main 
agronomic traits of D62B other than the expected BLB 
resistance. 

In contrast, there were evident morphological differ- 
ences between the two distinct rice varieties D62B and 
MH86. The sheath and blade of MH86 were light green 
but that of D62B were reddish dark green. Moreover, 
D62B showed a narrower blade, thinner and shorter 
stem, and a shorter growth period than MH86. 

Transcriptome profiling using Next-Gen sequencing 

In order to investigate transcriptome variations of the 
transgenic and MAB rice, digital gene expression profiles 
of DXT, DXB, D62B and MH86 were obtained using 
Next-Gen sequencing. The 5.8 to 6.1 million raw sequen- 
cing reads (Table 2) for the samples have been deposited 
into NCBI/GEO (accession number SRA061839). Prepro- 
cessing of these raw data (see Methods) gave rise to 5.7 to 
6.0 million qualified reads for down-stream analysis. The 
qualified reads were mapped to the cDNA sequences of 
Oryza Sativa Nipponbare with a stringent criterion, giving 
12.1 to 13.3 thousand unigenes per sample (Table 2). A 
statistical analysis (see Methods) of all the genes mapped 
in the four samples detected 15,268 unigenes, among 
which 10,362 were common to all four rice plants. Based 
on their abundance, measured in reads per kilobase of 
exon model per million mapped reads (RPKM), 508 
(4.9%), 3,786 (36.5%) and 6,068 (58.6%) of the 10,362 com- 
mon genes were expressed at high (RPKM >100), medium 
(10 < RPKM <100) and low (RPKM <10) levels, respect- 
ively. The results indicated that the majority of the 
expressed genes had low or median abundance. 

Transgenesis induced smaller transcriptome variation 
than MAB breeding 

A screening of differentially expressed genes (DEGs) re- 
vealed 984 DEGs (739 up-regulated and 245 down- 
regulated) in DXT vs. D62B, 1413 DEGs (938 up- and 
475 down-regulated) in DXB vs. D62B, and 2599 DEGs 
(2109 up- and 490 down-regulated) in MH86 vs. D62B 



Panicle length (cm) ^ Seed ratio (%) 1000-grains weight (g) 

23.9±1.2 94±3.0 25.8±2.8 

24.2±0.9 93±4.7 25.8±2.2 

23.8±1.1 93±2.9 26.0±1.7 

27.0±1.7* 80±5.8* 27.1 ±1.0 



Data included are averages of at least 5 Individual plants from which 3 infected leaves (a) or panicles (b) were scored. * represents a significant levels of p < 0.01 
beween D62B and MH86 (f-test). # represents the significant levels of p <0.01 in DXT and DXB with D62B as the control.D62B: the recipient control; DXT: 
transgenic line with Xall; DXB: MAB line with Xa21; IV1H86: another rice variety. 
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P1 P2 P3 P4 P6 P7 P8 P9 P10 P1 P2 P3 P4 P6 P7 P8 P9 P10 P1 P2 P3 P4 P6 P7 P8 P9 P10 




D62B 



DXT 



DXB 



Figure 1 Bacterial blighit resistance spectrum analysis of D62B, DXT and DXB. Nine Xoo strains from Pliilippine (Pl-PlO) were used to 
innoculate tlie rice leaves in tlie liigli tillering stage. The shorter the yellow portion on a leaf is, the higher plant resistance to Xoo infection. DXT: 
transgenic line with Xa21] DXB: MAB line with Xa21] D62B: the common parental line of DXT and DXB. 



(Additional file 1: Table SI). Eight up-regulated and two 
down-regulated genes in DXT were further analyzed by 
qRT-PCR (Additional file 2: Figure SI). Eight of these 
ten genes profiled exhibited expression consistent with 
the sequencing data, confirming the results from Next- 
Gen sequencing; seven of the eight up-regulated genes 
were indeed significantly up-regulated and one of the two 
down-regulated gene was significantly down-regulated in 
DXT (Additional file 2: Figure SI). The discrepancy be- 
tween the results of deep sequencing and qRT-PCR on 
two genes assayed may be due to the technical difference 
between the two techniques [24,25]. 

Remarkably, the transcriptome variations between the 
two rice varieties D62B and MH86, measured by the num- 
ber of DEGs between the two, was the largest, while the 
transcriptome variations between DXT and DXB was the 
smallest among the comparisons. More importantly, trans- 
genesis of Xa21 induced 43.6% less transcriptome variation 
than MAB breeding. The differences among such transcrip- 
tome variations were displayed in Figures 2A to 2C. Over- 
all, the results showed that the transcriptome change due 
to the introduction of Xa21 was significantly smaller than 
that by MAB breeding and the variation between the tran- 
scriptomes of two rice varieties D62B and MH86. 



Furthermore, a clustering analysis of DEGs (see Methods) 
showed that the rice plants with close genetic background 
tended to group together (Figure 2D). As expected, DXT 
and DXB were closer to their parental line D62B than they 
were to MH86. Despite the difference in the two breeding 
strategies, the transgenic line DXT was more closely re- 
lated to D62B than the MAB line DXB (Figure 2D). These 
results revealed that transgenesis and MAB breeding did 
not alter the transcriptome more significantly than an- 
other natural rice variety (i.e., MH86). As far as transcrip- 
tome variation is concerned, these results suggest that the 
transgenic rice is closer to its parental line than the MAB 
rice is to the same parental line. 

Transgenesis had less impact on molecular and cellular 
functions than MAB breeding 

In order to appreciate the possible consequences of intro- 
ducing Xa21 into D62B by the two breeding strategies, GO 
functional analysis was performed on the DEGs between 
DXT and D62B and between DXB and D62B (Figure 3). 
In total, 31 and 35 functional categories were enriched 
among the DEGs between DXT and D62B and between 
DXB and D62B, respectively. Surprisingly, 30 of these 
enriched functional categories were common (Figure 3B), 



Table 2 Statistics of RNA sequencing data, generated by lllumina sequencing, on four rice plants 



D62B 



DXT 



DXB 



IVIH86 



Raw reads 


5,801,661 


6,161,959 


6,101,257 


6,104,443 


Qualified reads 


5,677,504 


6,026,127 


5,962,165 


5,954,422 


Unambiguous tag mapping to Gene 


2,330,333 


2,724,696 


2,549,799 


2,665,554 


Unambiguous tag-mapped Genes 


12,048 


12,686 


12,526 


13,349 


Unambiguous tag/gene 


193.4 


214.8 


203.6 


199.7 



D62B: the recipient control; DXT: transgenic line with Xall; DXB: MAB line with Xall; MH86: another rice variety. 
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D62B 
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pXB 


O 


MH86 
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Figure 2 Transcriptome variations and expression relationships among four rice plants studied. Scatter plotes of transcriptome 
comparison of (A) DXT vs. D62B, (B) DXB vs. D62B, and (C) MH86 vs. D62B, where the horizontal and vertial axes represent the digital gene 
expression abundance after log 10 transformation, and dots in red represent differentially expressed genes. (D) Clustering results on the 
differentially expressed genes of DXT, DXB and MH86 with respect to D62B. Each column in the figure refers to a gene. Digital expression 
abundance, after log 10 transformation, is used in the plot. DXT: transgenic line with Xq21; DXB: MAB line with Xq21; D62B: the common parental 
line of DXT and DXB; MH86 :another rice variety. 



indicating that most of the DEGs in DXT and DXB have 
the same or similar functions. The single functional cat- 
egory specific to the DEGs in DXT is the "extracellular 
region part", and the four categories unique to the 
DEGs in DXB include "antioxidant", "enzyme regula- 
tor", "molecular transducer" and "death" (Figure 3B). 
Although it is unclear what phenotypic consequence, if 
any, that these DEGs may lead to, it is evident that the 
transcriptome variation caused by the DEGs in DXT 
imposed less cellular perturbation than that by the 
DEGs in DXB. In other words, the transgenesis of Xa21 
induced less disturbance to the molecular and cellular 
machineries than the MAB breeding did. 

To further understand the effects of the two different 
breeding strategies, we analyzed the genes expressed ex- 
clusively in DXT and exclusively in DXB, along with the 
ones expressed exclusively in D62B. We identified 758, 
821 and 404 genes that were exclusively expressed in 
DXT, DXB and D62B, respectively (Figure 4A). The 
DXT-, DXB- and D62B-specific genes were enriched in 
21, 33 and 26 GO functional categories, respectively 
(Figure 4B). Interestingly, the 21 categories enriched in 
the DXT-specific genes were also enriched in the DXB- 
and D62B-specific genes, and the 26 categories enriched 
in the D62B-specific genes were enriched in the DXB- 
specific genes as well (Figure 4B). This result indicated 



that although there were a large number of genes 
expressed exclusively in these plants, they carried the 
same or similar molecular and/or cellular functions, 
reflecting the functional elasticity and robustness of 
these genes. On the other hand, comparing with the 
genes expressed in the parental line, the genes expressed 
exclusively in DXT lost functions in 5 categories, while 
the genes expressed exclusively in DXB gained additional 
functions in 7 categories. It is possible that these lost 
and gained functions may be compensated for by the 
genes commonly expressed in these three rice plants. 
Nevertheless, this plant-specific gene function analysis 
also indicated that the perturbation to the molecular and 
cellular functions induced by transgenesis has a smaller 
or compatible scale than by MAB breeding. 

Genes exclusively expressed in transgenic rice functioned 
in pathogen defense 

Among the DXT-specific genes, 52 encoded transposons 
and retrotransposons, supporting the notion that mobile 
elements are typically activated during genetic trans- 
formation [26]. In addition, 40 DXT-specific genes 
encoded protein kinases, receptor-like protein kinases 
(RLKs), or OsWAK receptor-like protein kinases that 
have been implicated in stress/defense signaling-signal 
perception and signal transduction, which are consistent 
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Figure 3 Comparative analysis of the DEGs in DXT and DXB with respect to D62B. (A) Venn diagram for the numbers of DEGs in DXT and 
DXB. (B) GO functional categories of the DEGs in DXT and DXB. The y-axis on the right indicates the number of genes in a category, and the y- 
axis on the left is the percentage of genes to be analyzed in a category. DXT: transgenic line with Xq21] DXB: MAB line with Xq21] D62B: the 
common parental line of DXT and DXB. 



with the existing results that show that protein phos- 
phorylation and dephosphorylation play an important 
role in Xa21 induced gene regulation in response to 
pathogen invasion [27-32]. The remaining DXT-specific 
genes included 14 genes encoding transporter (including 
potassium, sulfate, metal cation, ABC, and ctr copper 
transportor), 16 genes encoding putative transcription 
factors (including 6 AP2 domain containing proteins, 1 
ethylene-responsive transcription factor (ERF), and 9 
MYB family transcription factors), 8 DXT-specific ex- 
pressed genes encoding resistance or resistance-like 
proteins (including 5 NBS-LRR type of proteins such as 
LRIO [33], RPMl [34] and I2GA-SH194-2 [35]) and 2 
genes encoding DNA methylation related proteins. DNA 
methylation has been reported to affect rice resistance 
response [36-38]. The AP2/ERF family is a large family of 
plant specific transcription factors that share a well- 
conserved DNA-binding domain that has been reported 
to activate the expression of abiotic stress-responsive 
genes via specific binding to the dehydration-responsive 



element/C-repeat (DRE/CRT) cis-acting element in their 
promoters [39]. The MYB family has key transcription fac- 
tors for controlling plant development and response to bi- 
otic and abiotic stresses [40] . Taken together, these results 
strongly suggested that transgenesis of Xa21 affected 
mainly those genes that were involved in pathways related 
to stress response. 

Compatible observations were also made on the genes 
expressed exclusively in DXB, i.e., these genes had the 
same or similar functions as those exclusively expressed 
in DXT. Among the DXB-specific genes, 58 encoded 
transposons and retrotransposons and 24 encoded pro- 
tein kinases, receptor-like protein kinases (RLKs) or 
OsWAK receptor-like protein kinases. The remaining 
DXB-specific genes included 10 genes that encode tran- 
sporters, 10 that encode putative transcription factors 
(including 3 MYB domain containing proteins), 13 that 
encode resistance or resistance-like proteins (including 7 
NBS-LRR type of proteins), and 3 that encode DNA 
methylation-related proteins. 
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Figure 4 Comparative analysis of the genes exclusively expressed in transgenic line DXT, MAB breeding line DXB and their parental 
line D62B. (A) Venn diagram for the numbers of genes exclusively expressed in DXT, DXB and D62B. (B) GO functional categories of the genes 
in DXT and DXB. The y-axis on the right indicates the number of genes in a category, and the y-axis on the left is the percentage of genes to be 
analyzed in a category. DXT: transgenic line with Xq21; DXB: MAB line with Xq21] D62B: the common parental line of DXT and DXB. 



Most pathways perturbed by transgenesis were also 
disturbed by MAB breeding 

Our final step of the comparative study of transcriptome 
variation was a pathways analysis on the DEGs in DXT 
and the DEGs in DXB. To this end, the DEGs in these 
two plants were mapped to the pathways in the Kyoto 
Encyclopaedia of Genes and Genomes (KEGG) database 
and searched for significantly enriched KEGG pathways 
that were potentially affected by gene expression variation 
[41]. Among the DEGs in DXT and DXB, 476 and 700 
were annotated on KEGG pathways, respectively. In 



comparison, more DEGs, i.e., 1331, between MH86 and 
D62B were mapped to KEGG pathways. With a statistical 
significance of j?<0.05; 11, 16, and 20 signaling pathways 
were abundant among the DEGs of DXT, DXB and 
MH86, respectively (Table 3, Additional file 1: Table S2). 
Seven of the 11 significantly enriched pathways in the 
transgenic rice were also significantly affected in the MAB 
rice. The rest four pathways (including Carotenoid biosyn- 
thesis. Regulation of autophagy, Arachidonic acid metab- 
olism and Anthocyanin biosynthesis) were also detected in 
the MAB rice but at a less significant level above the cutoff 



Gao et at. BMC Genomics 2013, 14:738 
http://www.bionnedcentral.conn/1471 -21 64/1 4/738 



Table 3 Pathways enriched in differentially expressed 
genes in the transgenic line DXT, MAB breeding line DXB 
and a rice variety MH86 with respect to recipient line 
D62B 



Pathways DXT DXB MH86 
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Metabolism of xenobiotics by 

r\/tnrhrnrnp P4Sn 
i^y n^^^i 1 1 i^ji lie tt^ju 




V 


Qnlirpn<;nrnp 

_)|-JIIL.C>^J>^l 1 IC 




V 


Endocytosis 




V 


Biosynthesis of plant hormones 




V 


Biosynthesis of unsaturated fatty acids 




V 


Glutathione metabolism 




V 


Linoleic acid metabolism 




V 


Biosynthesis of terpenoids and steroids 




V 


alpha-Linolenic acid metabolism 




V 


Phosphatidylinositol signaling system 




V 


Limonene and pinene degradation 




V 


Non-homologous end-joining 




V 


Biosynthesis of alkaloids derived from 
ornithine, lysine and nicotinic acid 




V 


Glycosylphosphatidylinositol 
(GPI)-anchor biosynthesis 




V 



Listed are significantly enriched KEGG terms with cutoff of p < 0.05. "V" means 
significantly affected pathway. Note that the four pathways enriched In DXT 
but not listed under DXB In the table were also enriched in DXB at a less 
significant level (with p-values from 0.06 to 0.88). 
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threshold (Additional file 1: Table S3). Specifically, "Carot- 
enoid biosynthesis", "Regulation of autophagy", "Arachi- 
donic acid metabolism", and "Anthocyanin biosynthesis" 
had /7-values of 0.22, 0.88, 0.06, and 0.67 respectively in 
the MAB rice. 

Discussion 

MAB breeding is a well accepted breeding technology 
for production of safe crops. Comparing crops produced 
from transgenesis and MAB breeding thus offers a viable 
approach to assessing the potential impact that trans- 
genic crops might have on human health. As a general 
principle, substantial equivalence [11] can be followed to 
assess the safety of transgenic crops. Using this principle 
as a guideline and yardstick, various aspects of substan- 
tial equivalence, e.g., molecular, biological, chemical and 
nutritional equivalence, can be examined. Taking a 
systems-biology perspective, we focused on molecular 
and biological equivalences between transgenic crops 
and plants from MAB breeding in the current study. 
This was done by contrasting transcriptome variations 
of a transgenic rice and a rice by MAB breeding. As a 
caveat, any results from such a study should not be 
interpreted to directly address the safety of transgenic 
plants. But rather, in addition to gaining insights, at a 
genome scale, into the biological processes and pathways 
that might be perturbed by transgenesis, such results 
can be used to pave the way for further study of chem- 
ical and nutritional equivalence. 

Adopting the above view point and taking Xa21, the 
most widely used resistance gene to rice bacterial leaf 
blight, as the target gene, we studied the phenotype and 
transcriptome differences between the transgenic and 
MAB rice carrying Xa21, Included in our expriment sys- 
tem were a transgenic rice (i.e., DXT) and a MAB rice 
(i.e., DXB) that were developed from the same parental 
line D62B. DXT was identical to D62B except the key 
Xa21 gene. The genetic makeup of DXB, through careful 
selection over six generations of backcrossing to D62B, 
was almost identical to that of D62B. Therefore, with 
exception of the specific breeding techniques used, DXT 
and DXB carried nearly identical genetic materials, mak- 
ing them ideal for the comparison of the two breeding 
techniques. A critical component of our experiments 
and analysis was another natural rice variety, MH86, 
which does not carry Xa21, Here, MH86 played the role 
of a"control", when compared with D62B, to reveal the 
natural range of transcriptome variation between the 
two natural rice plants. 

The comparison of transgenic and MAB rice was done 
at two levels. At the physiological level, we confirmed 
through extensive, although laborious, in-field screening 
of various morphological properties to ensure that DXT 
and DXB were morphologically equivalent or similar. 
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Due to genetic diversity, the morphological difference 
between natural varieties are typically greater than that 
of transgenesis, as we saw in our comparison between 
MH86 vs. D62B and DXT vs. D62B. At the molecular 
and cellular level, we carefully studied transcriptome 
variations as well as function and pathway alterations in 
DXT and DXB. The scale of transcriptome variations 
caused by transgenesis was smaller than that by MAB 
breeding. More importantly, the molecular and cellular 
functions that may be affected by transgenesis were also 
affected by MAB breeding, and the majority of pathways 
perturbed in the transgenic rice were also distorted in the 
MAB rice. In stark contrast, transcriptome and function 
variations between MH86 and D62B, which can be 
regarded as natural variations, were substantially larger 
than that between DXT and D62B. In short, the variations 
caused by transgenesis were smaller than that of MAB 
breeding and were within the natural range of variations. 

The results from the current study were in agreement 
with that in the literature on transgenic wheat and maize. 
For example, Gregersen et al compare the gene expression 
profiles of wildtype and transgenic wheat expressing an as- 
pergillus fumigatus phytase [42]. The results show that the 
expression profiles of the two plants are not significantly 
different. A similar study of transgenic wheat, which carries 
additional High Molecular Weight subunit genes, reaches 
the same conclusion that transgenesis has less impact 
on the transcriptome of wheat grain than conventional 
breeding [10]. Further, Coll et al uses microarrays to com- 
pare the expression profiles of commercial maize variety 
MON810 and near-isogenic varieties in leaves in vitro and 
also field cultured plants of AristisBt/Aristis and PR33P67/ 
PR33P66 [43,44]. The target gene is Bt, a gene non-native 
to maize encoding insecticidal crystal protein of the soil 
bacterium Bacillus thuringiensis (B.t), The results show 
that gene expression profiles of MON810 and comparable 
non-GM maize varieties are more similar to that of conven- 
tional lines and natural variation. 

On the other hand, the experiment designs of previous 
studies are far from perfect. In [10], comparative analysis 
of the transcriptomes of three wheat lines, including a 
parental line, a conventional breeding line and a trans- 
genic line, reaches the conclusion that transgenesis has 
less impact on the transcriptome of wheat grain than 
conventional breeding. However, several additional genes 
or sequences, including the marker gene (Bar gene), the 
reporter gene {uidA gene), and sequences derived from 
the bacterial plasmid, are present in the transgenic line. 
The transgenic line, referred to as a "clean fragment" 
line, still contains the bar gene. In fact, selective 
markers, such as ban are the main concern of biosafety 
of transgenic crops which induce specific pleiotropic ef- 
fects [45]. In the studies reported in [42-44], foreign 
transgenes are used. The current study overcame these 



issues. First, the target gene Xa21 was cloned from 
Oryza longistaminata, which is native to rice. Second, 
other than Xa21 sequence, DXT did not carry any add- 
itional sequences such as that for marker genes and re- 
porter genes. Third, unlike the previous studies where a 
target gene is under the control of a foreign promoter 
[42,46], in the current study Xa21 in DXT carried its 
own promotor. As a result of these factors, the results 
from our analysis were specific to the functions of Xa21 
rather than some artifacts of genes foreign to rice. 

More importantly, our analysis expanded beyond tran- 
scriptome profiling to include analyses of gene functions 
and signaling pathways that might be altered by the 
introduction of Xa21, Our results showed that most of 
the molecular and cellular functions affected by trans- 
genesis were influenced by MAB and functional categor- 
ies that affected by MAB were more than those that 
affected by transgenesis. Analyzing the pathways per- 
turbed by transgenesis and MAB showed that majority 
pathways altered in transgenic rice were also distorted in 
MAB rice. The bigger difference between the transcrip- 
tome variations of DXB and D62B could be attributed to 
two factors. Firstly, the transgenic rice analyzed in our 
study was carefully selected to have favorable properties 
for the purpose of the study. In particular, DXT was a 
single-copy, marker-free, and homozygous Xa21 trans- 
genic line and also had consistent agronomic traits simi- 
lar to that of the parental line D62B. The genetic 
makeup of DXT was identical to D62B with the excep- 
tion of the Xa21 gene. Secondly, although six back- 
crosses have been used to eliminate linkage drag, a few 
unwanted donor chromosomal segments could still have 
been retained in DXB when Xa21 was introgressed into 
D62B. Therefore, the transcriptome of DXT was closer 
to that of D62B than DXB to D62B. 

Thanks to Next-Gen sequencing, the gene profiling 
experiments of our study were genome wide and pro- 
vided a high resolution [47]. In contrast, the previous 
studies focus on a limited number of genes using micro- 
array profiling which is less accurate and restricted to 
annotated genes [10,42-44,46,48]. Such a large scale 
functional analysis performed in the current study has 
never been attempted before and the functional analyses 
provided deep insight into the functional equivalence be- 
tween the transgenesis and MAB breeding. 

Conclusions 

MAB breeding and transgenesis are two most popular 
breeding techniques for producing plants of favorable 
traits. As a newer biotechnology, transgenesis has made 
visible contributions to increase yield of staple crops. 
While plants produced by MAB breeding have already 
been widely accepted, transgenic plants are facing the 
challenge regarding safety on human health. In the 
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current study, these two distinct breeding techniques 
was closely compared. The analysis focused on tran- 
scriptome variations in rice plants generated from these 
two breeding techniques and on that of two natural rice 
plants as a baseline of the comparison. The study com- 
bined careful assessment of agronomic traits, transcrip- 
tome profiling by Next-Gen sequencing, and functional 
and pathway analyses. Two important conclusions can 
be drawn from the results. First, transcriptome variation 
caused by transgenesis is significantly smaller than that 
by MAB breeding and is within the range of natural vari- 
ation. Second, the functional categories of differentially 
expressed genes due to these two breeding techniques and 
the pathways perturbed by these techniques are not sub- 
stantially different. These results suggest the transgenic 
rice and rice from MAB breeding that were compared in 
the current study are substantially equivalent at the mo- 
lecular and biological levels. The data and results can be 
used to study chemical and nutritional equivalence of rice 
generated by transgenesis and MAB breeding. 

Methods 

Resistance analysis of rice varieties 

Both D62B and MH86 are indica rice and widely used 
as parental lines of hybrid rice in China. The transgenic 
line DXT with Xa21 was developed from D62B in our 
previous study [22]. The MAB line DXB with Xa21 was 
developed by backcrossing to D62B in the current study. 
When the rice grew to the high tillering stage, five to 
seven fully expanded leaves were inoculated with Xoo 
using the leaf-clipping method [49]. The cultures were 
grown on PSA (Potato-Sugar- Agar) medium (potato, 300 
g/L; Ca(N03)2*4H20, 0.5 g/L; Na2HP04^12H20, 2.0 g/ 
L; sugar, 15 g/L; agar 15 g/L) at 28°C for 3 days. Inocu- 
lums were prepared by suspending the bacterial cells in 
sterile water and adjusting the concentration to about 
10^ cells per milliliter. Phenotype scoring was carried 
out at 15 days post innoculation (dpi). 

Molecular analysis of transgenic and MAB 
breeding plants 

Genomic DNA was isolated from fresh rice leaf tissue 
using the cetyltrimethylammonium bromide protocol. 
Transgenic plants were validated by polymerase chain 
reaction (PGR) and Southern blotting. The primers 
which were used for molecular analysis of transgenic 
plants were described detailedly in Gao et al. (2011). 
MAB breeding plants were validated by PGR using the 
same primers as transgenic plants. 

Quantitative PGR 

DNasel-treated RNA was used for fist strand cDNA syn- 
thesis using M-MLV reverse transcriptase (Promega) and 
oligo(dT)i5 according to the manufactures protocols. 



Specific pairs of primers for SYBR-green detection and 
quantification of selected DEGs were designed using Pri- 
mer Express software follwing the primer design guidance. 
The primer sequences are listed in Additional file 1: Table 
S5. The endogenous actin gene was run in parallel as con- 
trol PGR reaction and the untransfromed receptor line was 
used as the calibrator to normlize the relative expression 
levels of the target DEG. Triplicate samples for each tested 
line were prepared for real-time PGR assays. The 2'^^^^ 
method was used to calculate relative changes of gene ex- 
pression [50]. PGR reactions were performed using the 
TranStart Green qPGR SupMix reagent (TransGen Biotech, 
Inc.) on BioRad GFX96 PGR system. 

Illumina sequencing 

Sequencing and preliminary data acquisition were fin- 
ished by Beijing Genomics Institute. The experimental 
process includes sample preparation and sequencing. 
The main reagents and supplies used included Illumina 
Gene Expression Sample Prep Kit and Solexa Sequen- 
cing Ghip (flowcell), and the main instruments used in- 
cluded Illumina Gluster Station and Illumina HiSeq™ 
2000 Systerm. Ten leaves randomly selected from ten in- 
dividual plants of each tested line were harvested and 
pooled together for RNA extraction. Total RNA was iso- 
lated using Trizol reagent (Invitrogen Gorporation, 
Garlsbad, GA, USA) following the manufactures protocols 
and purified using an RNeasy Plant Kit (Qiagen, Hilden, 
Germany). The integrity and purity of RNA samples were 
determined by gel electrophoresis and OD 260/280 nm 
absorption ratios and the RNA concentration was quanti- 
fied using a NanoDrop ND-1000 spectrophotometer 
(NanoDrop Technology, Wilmington, DE, USA). 

Processing of sequencing data 

Raw sequencing data in image were transformed by base 
calling into raw reads. Raw reads were then transformed 
into clean reads after removal of such reads as 3' adaptor 
sequences, empty reads, low quality reads, reads which 
were too long or too short, and reads with only one copy 
(probably due to sequencing error). Glean reads that 
were mapped to reference sequences from multiple 
genes were removed. The remaining clean reads were 
designated as unambiguous clean or qualified reads. 

Gene expression and differentially expressed genes 

The number of qualified reads for each gene was tallied 
and normalized to TPM (number of transcripts per million 
qualified reads), which was then used as the digital gene ex- 
pression abundance of the gene. Genes that were differen- 
tially expressed, referred to as DEGs, across two plants 
were identified using the criterion of at least 2-fold change 
to digital expression abundance and with a cut-off of False 
Discovery Rate (FDR) no greater than 0.001. 
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Cluster analysis 

All DEGs in DXT, DXB and MH86, with respect to 
D62B, were aggregated. The logio-transformed TPM 
values of the genes in the combined set of DEGs were 
used to cluster the DEGs by hierarchical clustering with 
euclidean distance and average linkage. The clustering 
was done and the result was plotted using Multi Experi- 
ment Viewer (MeV) [51]. 

Gene ontology and KEGG pathway analyses 

GO functional category analysis was performed separ- 
ately on DEGs and genes exclusively expressed in a spe- 
cific plant, e.g., DXT or DXB. This analysis was done by 
mapping the set of genes of interest (DEGs or plant- 
specific genes) to the terms in plant GOSlim Ontologies 
from the Rice Genome Annotation Project (http://rice. 
plantbiology.msu.edu/annotation_pseudo_goslim.shtml). 
WEGO (Web Gene Ontology Annotation Plot [52]) tool 
was used to plot GO annotation results. For pathway en- 
richment analysis, all DEGs involved in a comparison 
were mapped to the terms in the KEGG database to 
identify significantly enriched KEGG terms. The Kobast 
2.0 [53] tool for pathway enrichment analysis was used 
with a cutoff of p value no greater than 0.05. 

Availability of supporting data 

The gene expression profiling data on the four rice 
plants from Illumina HiSeq™ 2000 has been deposited 
into NCBI/GEO database under the accession number 
SRA061839. 
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